home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / map_set.pro < prev    next >
Text File  |  1997-07-08  |  39KB  |  981 lines
  1. ; $Id: map_set.pro,v 1.64 1997/04/08 21:29:26 dave Exp $
  2. ;
  3. ; Copyright (c) 1993-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5.  
  6. PRO map_struct_append, In_struct, Name, Value, SUPERCEDE = super
  7. ;
  8. ; Append/Replace the tagname Name, with a given value to the In_Struct.
  9. ; If SUPERCEDE is set, and a tag with the given Name already exists, 
  10. ;    replace its value.
  11. ; If SUPERCEDE is NOT set, and the given tag exists, do nothing.
  12. ;
  13.  
  14. ntags=N_TAGS(In_struct)
  15.  
  16. if ntags eq 0 then begin        ;If In_struct is undef, just return the tag
  17.     In_struct = CREATE_STRUCT(name, value)
  18. endif else begin
  19.     tnames=tag_names(IN_STRUCT)
  20.     index=where(tnames eq NAME, count)
  21.     if count eq 0 then $        ;No match, append
  22.       In_Struct = CREATE_STRUCT(in_struct, name, value)    $
  23.     else if keyword_set(super) then $
  24.       in_struct.(index[0]) = value ;Overwrite value?
  25.                                 ;Otherwise, tag is already there, don't add
  26. endelse
  27. end
  28.  
  29.  
  30.  
  31. PRO MAP_STRUCT_MERGE, orig, add, SUPERCEDE = super
  32. ; Add the keywords in the struct add to the original struct orig.
  33. ; ADD must be either undefined or defined as a structure.
  34. ; If a given tag name is present in both structures:
  35. ;   If SUPERCEDE is set, then the value of add.tag replaces the
  36. ;   original value.  If SUPERCEDE is not set, the original value
  37. ;   remains unchanged.
  38.  
  39. if n_elements(orig) eq 0 then begin
  40.     if n_elements(add) ne 0 then orig = add
  41.     return
  42. endif
  43.  
  44. otag = tag_names(orig)
  45. atag = tag_names(add)
  46. for i=0, n_elements(atag)-1 do begin ;Add the new ones
  47.     match = where(otag eq atag[i], count)
  48.     if count eq 0 then $
  49.       orig = create_struct(orig, atag[i], add.(i)) $ ;Disjoint tag name, add it
  50.     else if keyword_set(super) then $ ;Got a match, if super is set, overwrite.
  51.       orig.(match[0]) = add.(i) ;Common, add if set
  52. endfor
  53. end
  54.  
  55.  
  56. FUNCTION map_point_valid, lon, lat, u, v
  57. ; Return 1 if the point at geo-coordinates (lon, lat) is mappable,
  58. ;    0 if not.
  59. ; Mappable means that the point is transformable, and is not clipped
  60. ;   by the mapping pipeline.
  61. ; Input: lon, lat = input coordinates in degrees.
  62. ; Output: u, v = uv coordinates if point is mappable.  The contents of
  63. ;   u and v are undefined if the function returns 0, indicating that the
  64. ;   point is not mappable.
  65. ; 
  66.  
  67. u = 0.0 & v = 0.0
  68. p = convert_coord(lon, lat, /DATA, /TO_NORM) ;Cvt to normalized coords
  69. if finite(p(0)) eq 0 then return, 0
  70. u = (p(0) - !x.s(0)) / !x.s(1)  ;To UV coords
  71. v = (p(1) - !y.s(0)) / !y.s(1)
  72.  
  73. clat = cos(lat * !dtor)         ;Cvt lon lat to cartesian 3D coords
  74. xyz = [ clat * cos(lon * !dtor), clat * sin(lon * !dtor), sin(lat * !dtor)]
  75.  
  76. p = !map.pipeline
  77. s = size(p)
  78. for i=0, s(2)-1 do begin        ;# of stages
  79.     case p(0,i) of
  80.         0 : return, 1           ;0 = END, means we're successful.
  81.         1 : dummy = 0           ;Ignore splits
  82.         2 : if total(xyz * p(1:3,i)) + p(4,i) lt 0 then $ ;Outside clip plane?
  83.           return, 0
  84.         3 : dummy = 0           ;Ignore transforms
  85.         4 : if u lt p(1,i) or v lt p(2,i) or $ ;Within UV bounding box?
  86.           u gt p(3,i) or v gt p(4,i) then return, 0
  87.     endcase
  88. endfor
  89. return, 1                       ;Should never hit here.
  90. end                             ;map_point_valid
  91.  
  92.  
  93.  
  94. PRO map_set_ll_box
  95.  
  96. ;Try to determine a lon/lat range from the existing parameters, and
  97. ;save in !MAP.LL_BOX = [latmin, lonmin, latmax, lonmax].
  98. ; If we can't reliably determine the lat/lon limits then the
  99. ; coorresponding min and max are both set to zero.
  100.  
  101. ; If either pole is mappable, then all longitudes are mappable
  102. north = map_point_valid(0, 90)
  103. south = map_point_valid(0, -90)
  104. dateline = (north + south eq 0) and map_point_valid(-180., !map.p0lat) and $
  105.   (map_point_valid(0., !map.p0lat) eq 0) ; Dateline visible
  106.  
  107. if north and south then return  ;The whole globe is visible
  108.  
  109. if north then begin
  110.     lonmax = 180. & lonmin = -180. ;All longitudes are visible
  111.     latmax = 90. & latmin = 90.
  112. endif else if south then begin
  113.     lonmax = 180. & lonmin = -180.
  114.     latmax = -90. & latmin = -90.
  115. endif else begin
  116.     if dateline then begin
  117.         lonmax = 0. & lonmin = 360.
  118.     endif else begin
  119.         lonmax = -180. & lonmin = 180.
  120.     endelse
  121.     latmin = 90. & latmax = -90.
  122. endelse
  123.  
  124. box = !map.uv_box
  125.  
  126. ; Edges are in order: bottom, right, top, left.
  127.  
  128. sidex = [0, 2, 2, 0 ,0]         ;X index of UV corner
  129. sidey = [1, 1, 3, 3, 1]         ;Y index
  130. uvcen = (box([0,1]) + box([2,3]))/2.
  131. del = 8                         ;# of intervals
  132. found = bytarr(4)
  133.  
  134. ; March along the 4 uv edges of the plot, keeping the range in
  135. ; geo-coordinates.
  136.  
  137. for iside = 0, 3 do begin
  138.     s0 = box([sidex(iside), sidey(iside)]) ;Start of edge
  139.     s1 = box([sidex(iside+1), sidey(iside+1)]) ;End of edge
  140.     for i = 0, del-1 do begin
  141.         s2 = s0 + (s1-s0) * (float(i)/del) ;UV coordinate
  142.         n2 = s2 * [!x.s(1), !y.s(1)] + [!x.s(0), !y.s(0)] ;Normalized coord
  143.         p = convert_coord(n2(0), n2(1), /NORMAL, /TO_DATA)
  144.         if finite(p(0)) then begin
  145.             if dateline and (p(0) lt 0) then p(0) = p(0) + 360
  146.             lonmax = lonmax > p(0)
  147.             lonmin = lonmin < p(0)
  148.             latmax = latmax > p(1)
  149.             latmin = latmin < p(1)
  150.             found(iside) = 1    ;Found a legit point
  151.         endif
  152.     endfor                      ;for i
  153. endfor                          ;for iside
  154.  
  155. if total(found) eq 4 then begin ;If we didn't values along 4 sides then punt
  156.     !map.ll_box = [latmin, lonmin, latmax, lonmax]
  157. ;    print, !map.ll_box, format='(4f10.2)'
  158. endif
  159. end
  160.  
  161.  
  162.  
  163. PRO map_set_limits, limit, uvrange
  164. ; Limit = 4 or 8 point vector containing lat/lon limits. (input).
  165. ;    WARNING: limit may be changed.
  166. ; uvrange = uv extent of map in form of [umin, vmin, umax, vmax]
  167. ; (output).  Also, on output, limit is converted to the 8 point form.
  168. ;
  169.  
  170. ; Set unit scaling for convert_coord:
  171. !x.s = [0,1]
  172. !y.s = [0,1]
  173.  
  174. if n_elements(limit) eq 4 then begin ;cvt 4 pt to 8 pt limit.
  175. ; ***********        limit = [latmin, lonmin, latmax, lonmax]
  176.     n = 7                       ;Check an NxN grid within the lat/lon limits
  177.     londel = limit[3] - limit[1]
  178.     if londel lt 0 then londel = londel + 360
  179.     lats = findgen(n) * ((limit[2]-limit[0])/(n-1)) + limit[0]
  180.     lons = findgen(n) * (londel / (n-1)) + limit[1]
  181.     j = 0
  182.     for i=0, n^2-1 do begin     ;Check the grid for extrema
  183.         uv = convert_coord(lons[i mod n], lats[i/n], /TO_NORM, /DATA)
  184.         if finite(uv[0]) then begin
  185.             if j ne 0 then begin ;Already got a limit?
  186.                 uvrange[0] = uvrange[0] < uv[0]
  187.                 uvrange[2] = uvrange[2] > uv[0]
  188.                 uvrange[1] = uvrange[1] < uv[1]
  189.                 uvrange[3] = uvrange[3] > uv[1]
  190.             endif else uvrange = uv[[0,1,0,1]]
  191.             j = j + 1
  192.         endif                   ;in range
  193.     endfor
  194.     if j lt 2 then message, 'At least two limit points must be mappable', /INFO
  195. endif else if n_elements(limit) eq 8 then begin
  196. ; 8 point limit in form of:
  197. ;    [latLeft,lonLeft, latTop, lonTop, LatRt, lonRt, LatBot, LonBot]
  198. ;
  199.     for i=0,3 do begin
  200.         uv = convert_coord(limit[i*2+1], limit[i*2], /TO_NORM, /DATA)
  201.         if finite(uv[0]) eq 0 then $
  202.           message, 'Unmappable limit point: ' + $
  203.             string(limit[i*2+1], limit[i*2]), /INFO
  204. ;    Left = uvrange(0)=uv(0); Top = uvrange(3) = uv(1); Right =
  205. ;    uvrange(2) = uv(0); Bot = uvrange(1) = uv(1)
  206.         uvrange[([0,3,2,1])[i]] = uv[([0,1,0,1])[i]]
  207.     endfor
  208. endif else message, 'Map limit must have 4 or 8 points'
  209. end
  210.  
  211.  
  212.  
  213. FUNCTION map_rotxyz, p, rx, ry, rz ;Rotate the vector p(3) counterclockwise
  214. ; about the X, Y, and Z, by the amounts rx, ry, rz, degrees, in order.
  215. ;
  216. p1 = p
  217. dtor = !dpi/ 180.
  218.  
  219. if rx ne 0.0 then begin
  220.     sx = sin(rx * dtor) & cx = cos(rx * dtor)
  221.     t = [[1,0,0],[0,cx, sx], [0, -sx, cx]]
  222.     p1 = t # p1
  223. endif
  224. if ry ne 0.0 then begin
  225.     sy = sin(ry * dtor) & cy = cos(ry * dtor)
  226.     t = [[cy, 0, -sy], [0, 1, 0], [sy, 0, cy]]
  227.     p1 = t # p1
  228. endif
  229. if rz ne 0.0 then begin
  230.     sz = sin(rz * dtor) & cz = cos(rz * dtor)
  231.     t = [[cz, sz, 0], [ -sz, cz, 0], [0,0,1]]
  232.     p1 = t # p1
  233. endif
  234. return, p1
  235. end
  236.  
  237.  
  238. FUNCTION GREAT_CIRCLE, lon0, lat0, lon1, lat1, RADIANS=radians, PRINT=print
  239. ; Return: [sin(c), cos(c), sin(az), cos(az)] from the great circle
  240. ; distance between two points.  c is the great circle angular distance,
  241. ; and az is the azimuth between two points.
  242.  
  243. if KEYWORD_SET(radians) eq 0 then k = !dtor else k = 1.0
  244.  
  245. coslt1 = cos(k*lat1)
  246. sinlt1 = sin(k*lat1)
  247. coslt0 = cos(k*lat0)
  248. sinlt0 = sin(k*lat0)
  249.  
  250. cosl0l1 = cos(k*(lon1-lon0))
  251. sinl0l1 = sin(k*(lon1-lon0))
  252.  
  253. cosc = sinlt0 * sinlt1 + coslt0 * coslt1 * cosl0l1 ;Cos of angle between pnts
  254. sinc = sqrt(1.0 - cosc^2)
  255. if abs(sinc) gt 1.0e-6 then begin 
  256.     cosaz = (coslt0 * sinlt1 - sinlt0*coslt1*cosl0l1) / sinc ;Azmuith
  257.     sinaz = sinl0l1*coslt1/sinc
  258. endif else begin        ;Its antipodal
  259.     cosaz = 1.0
  260.     sinaz = 0.0
  261. endelse
  262.  
  263. if keyword_set(print) then begin
  264.     print, 'Great circle distance: ', acos(cosc) / k
  265.     print, 'Azimuth: ', atan(sinaz, cosaz)/k
  266. ;    print,'sinaz, cosaz = ', sinaz, cosaz
  267. endif
  268. return, [sinc, cosc, sinaz, cosaz]
  269. end
  270.  
  271.  
  272. PRO map_satellite_limit, n, xr, yr
  273. ; Return n or less points, defining the limit of the limb of a satellite
  274. ; projection
  275. a = findgen(n) * (2 * !pi / (n-1))
  276. st = !map.p[0] > 1.0
  277. st = sqrt((st-1.)/(st+1.))    ;Map limit radius for satellite projection
  278. xr = cos(a) * st        ;Taken from Snyder Pg. 175
  279. yr = sin(a) * st
  280. a =  yr * !map.cosr + xr * !map.sinr
  281. xr = xr * !map.cosr - yr * !map.sinr
  282. yr = a
  283. a = yr * !map.p[2]/(!map.p[0]-1.0) + !map.p[3]
  284. good = where(a gt 0.0, count)
  285. if count eq 0 then begin        ;If nothing, return nothing
  286.     trash = temporary(xr) & trash = temporary(yr)
  287.     return
  288. endif
  289. if count lt n then begin        ;Extract & shift so pnts on map are contiguous
  290.     i0 = 0            ;Get 1st good pnt following a bad pnt
  291.     for i=1,count-1 do if good[i-1]+1 ne good[i] then i0 = i
  292.     good = shift(good,-i0)
  293. endif
  294. a = a[good]
  295. xr = xr[good] * !map.p[3] / a
  296. yr = yr[good]/a
  297. end
  298.  
  299.  
  300.  
  301. PRO map_proj_info, iproj, CURRENT=current, $
  302.        UVRANGE=r, WIDTH=meters, NAME=name, CYLINDRICAL=cyl, $
  303.        PROJ_NAMES = proj_names_p, CIRCLE=is_circular, $
  304.        AZIMUTHAL=is_azimuthal, UV_LIMITS = uv_limits, $
  305.        LL_LIMITS = ll_limits
  306. ; Return information about maps in general, or projection iproj.
  307.  
  308. ; Iproj = projection number, Input.  If CURRENT is set, then the
  309. ;     current map projection is returned in Iproj.
  310. ; CURRENT = if set, use current proj number, and return that index in iproj.
  311. ; UVRANGE =  UV limits of the given projection, [umin, vmin, umax, vmax],
  312. ; WIDTH = (output) the width, in meters, of the the map at a scale
  313. ;  of 1:1.  The radius of the earth is used. If this result is 0, the
  314. ;  projection's UV coordinates are already in meters.
  315. ; NAME = the name of the projection.
  316. ; CYLINDRICAL = 1 if the projection is cylindrical or pseudo-cylindrical.
  317. ; PROJ_NAMES = string array containing the names of all the
  318. ;     available projections.
  319. ; UV_LIMITS = returned rectangle of the current map in UV coordinates,
  320. ;    [Umin, Vmin, Umax, Vmax].
  321. ; LL_LIMITS = returned rectangle of the current map in latitude/longitude
  322. ;     coordinates, [Lonmin, Latmin, Lonmax, Latmax].  Note, the
  323. ;     range may not always be available, and the min and max may
  324. ;     both be set to 0.
  325.  
  326. ; Issues for a new projection:
  327. ; UVrange, proj_scale, is_cyl, is_azimuthal, circle, name, horizon, clipping
  328.  
  329. ; 
  330. common map_set_common, proj_scale, proj_is_cyl, proj_names, circle, $
  331.   proj_is_azim, user_defined
  332.  
  333. on_error, 2
  334. if n_elements(proj_scale) eq 0 then begin ;Init constant arrays?
  335.     p2 = 2 * !pi
  336.     s8 = 2 * sqrt(8)
  337.     ro = 2 * 2.628
  338.     
  339. ;                St Or LC La Gn  Az Sa  Cy  Me  Mo  Si  Ai  Ha Al Tm  Mi Ro
  340. ;                    1  2  3  4  5   6  7   8   9  10  11  12  13 14 15  16 17
  341.   proj_scale =  [ 0, 4, 2, 4, 4, 4, p2, 1, p2, p2, s8, p2, p2, s8, 4, 0, p2,ro]
  342.   proj_is_cyl = [ 0, 0, 0, 0, 0, 0,  0, 0,  1,  1,  1,  1,  1,  1, 0, 1,  1, 1]
  343.   proj_is_azim= [ 0, 1, 1, 0, 1, 1,  1, 0,  0,  0,  0,  0,  0,  0, 0, 0,  0, 0]
  344.   circle =      [ 0, 1, 1, 0, 1, 1,  1, 1,  0,  0,  1,  0,  1,  1, 0, 0, 0,  0]
  345.   proj_names =  $
  346.     ["", "Stereographic", "Orthographic", "LambertConic", "LambertAzimuthal", $
  347.      "Gnomic", "AzimuthalEquidistant", "Satellite", "Cylindrical", $
  348.      "Mercator", "Mollweide", "Sinusoidal", "Aitoff", "HammerAitoff", $
  349.      "AlbersEqualAreaConic", "TransverseMercator", $
  350.      "MillerCylindrical", "Robinson"]
  351.   user_defined = bytarr(n_elements(proj_names))
  352. endif
  353.  
  354. max_proj = n_elements(proj_names)-1
  355. if arg_present(proj_names_p) then proj_names_p = proj_names
  356. if KEYWORD_SET(current) then iproj = !map.projection ;Get current projection
  357.  
  358. if ARG_PRESENT(uv_limits) then uv_limits = !map.uv_box
  359. if ARG_PRESENT(ll_limits) then ll_limits = !map.ll_box
  360.  
  361. if n_elements(iproj) eq 0 then return ;Rest of params are proj specific
  362.  
  363. if iproj lt 1 or iproj gt max_proj then $
  364.   message, 'Projection number must be within range of 1 to'+ strtrim(max_proj,2)
  365.  
  366. ones = [-1., -1., 1., 1.]
  367. r_earth = 6378206.4d0 ;Earth equatorial radius, meters, ; Clarke 1866 ellipsoid
  368.  
  369. cyl = proj_is_cyl[iproj]
  370. name = proj_names[iproj]
  371. meters = proj_scale[iproj] * r_earth
  372. is_circular = circle[iproj]
  373.  
  374. if arg_present(r) then case iproj of
  375.     1: r = ones                 ;Stereo
  376.     2: r = ones                 ;Ortho
  377.     3: r = 2 * ones             ;Lambert Conic
  378.     4: r = 2 * ones             ;lambert
  379.     5: r = 2 * ones             ;gnomic
  380.     6: r =  !DPI * ones         ;azimuth
  381.     7: if !map.p[1] ne 0 then begin ;Complicated satellite?
  382.         map_satellite_limit, 60, xr, yr
  383.         if n_elements(xr) eq 0 then r = [0,0,0,0] $
  384.         else r = [min(xr, max=maxxr), min(yr, max=maxyr), maxxr, maxyr]
  385.         meters = (r[2]-r[0]) * r_earth
  386.     endif else begin
  387.         st = !map.p[0] > 1.0   ;For satellite
  388.         r = sqrt((st-1)/(st+1)) * ones ;Simple case satellite
  389.         meters = (r[2]-r[0]) * r_earth
  390.     endelse
  391.     8: r = [-1., -.5, 1., .5] * !DPI ;Cylind equidistant
  392.     9: r = [-!DPI, -2.43, !DPI, 2.43] ;Mercator  (+/- 80 degs of lat) ~ 98%
  393.     10: r = sqrt(2.0) * [-2.,-1.,2.,1.] ;Mollweide
  394.     11: r = [-1., -.5, 1., .5] * !DPI ;sinusoidal
  395.     12: r = [-1., -.5, 1., .5] * !DPI ;Aitoff
  396.     13: r = sqrt(2.0) * [-2., -1., 2., 1.] ; Hammer-aitoff
  397.     14: r = 2 * ones            ;Albers Conic
  398.     15: r = 2*!pi * r_earth * ones ;Transverse Mercator, Units = meters.
  399.     16: r = [-!DPI, -1.832, !DPI, 1.832] ;Miller Cylind  (+/- 80 degs of lat)
  400.     17: r = [-2.628, -1.349, 2.628, 1.349] ;Robinson
  401. endcase
  402.  
  403. end
  404.  
  405.  
  406. PRO map_horizon, FILL=fill, NVERTS=n, ZVALUE=zvalue, _EXTRA=Extra
  407. ; Draw/fill the horizon (limb) of a map
  408.  
  409. map_proj_info, /CURRENT, NAME=name, UVRANGE=p, CIRCLE=is_circular
  410.  
  411. r = !map.rotation        ;Our rotation
  412.                                 ;Can't do conics
  413. if name eq "LambertConic" or name eq "AlbersEqualAreaConic" then return
  414. if n_elements(n) le 0 then n = 60 ;# of vertices
  415. a = findgen(n+1) * (2 * !pi / n) ;N angles from 0 to 2 pi
  416.  
  417. if name eq "Satellite" and !map.p[1] ne 0 then begin ;Satellite proj, tilted
  418.     map_satellite_limit, 60, xr, yr
  419.     r = 0.0                     ;No more rotation
  420. endif else if is_circular then begin ;Std circular projection
  421.     xr = (p[2]-p[0])/2. * cos(a)
  422.     yr = (p[3]-p[1])/2. * sin(a)
  423. endif else if name eq "Sinusoidal" then begin ;Sinusoidal
  424. ;    flon = (!map.out(3)-!map.out(2))/360.
  425.     flon = 1.0
  426.     xr = flon * !pi * cos(a)
  427.     yr = fltarr(n+1)
  428.     for i=0,n do case fix(a[i]/(!pi/2)) of
  429.         0:    yr[i] = a[i]
  430.         1:    yr[i] = !pi-a[i]
  431.         2:    yr[i] = !pi-a[i]
  432.         3:    yr[i] = a[i]-2*!pi
  433.         4:    yr[i] = 0.0     ;Last pnt
  434.     endcase
  435. endif else if name eq "Robinson" then begin ;Robinson
  436.     lp = [2.628, 2.625, 2.616, 2.602, 2.582, 2.557, 2.532, 2.478, 2.422, $
  437.           2.356, 2.281, 2.195, 2.099, 1.997, 1.889, 1.769, 1.633, 1.504, 1.399]
  438.     lm = [0., 0.084, 0.167, 0.251, 0.334, 0.418, 0.502, 0.586, 0.669, 0.752, $
  439.           0.833, 0.913, 0.991, 1.066, 1.138, 1.206, 1.267, 1.317, 1.349]
  440.     xr = [lp, -reverse(lp), -lp, reverse(lp)]
  441.     yr = [lm, reverse(lm), -lm, -reverse(lm)]
  442. endif else if name eq "UserDefined" then begin ;User projection
  443.  
  444. endif else begin                ;Rectangular
  445.     xr = [p[0], p[2], p[2], p[0], p[0]]
  446.     yr = [p[1], p[1], p[3], p[3], p[1]]
  447. endelse
  448.  
  449. if r ne 0.0 then begin
  450.     t1 = xr * !map.cosr + yr * !map.sinr ;Now rotate by - angle...
  451.     yr =  yr * !map.cosr - xr * !map.sinr
  452.     xr = t1
  453. endif
  454.  
  455. xtsave = !x.type
  456. if n_elements(zvalue) eq 0 THEN zvalue = 0
  457.  
  458. !x.type=0            ;Plot in UV space
  459. if keyword_set(fill) then polyfill, xr, yr, _EXTRA=Extra, NOCLIP=0 $
  460. else plots, xr, yr, zvalue, _EXTRA=Extra, NOCLIP=0
  461. !x.type=xtsave
  462. end
  463.  
  464.  
  465. PRO map_set_split, ADD=add, EQUATOR = equator
  466. ;  Set the splitting hemisphere.  If EQUATOR is set, set the splitting 
  467. ;  point as if the center of projection was on the equator, i.e. use
  468. ;  a latitude of 0.
  469.  
  470. pole = !map.pole[4:6]        ;Locn of pole
  471.  
  472. if keyword_set(equator) then begin ;use lat = 0
  473.     split = [!map.p0lon, 0, -sin(!map.u0), cos(!map.u0), 0.0, 0.0]
  474. endif else begin
  475.     coslat = cos(!map.v0)
  476.     p0 = [ cos(!map.u0) * coslat, sin(!map.u0) * coslat, sin(!map.v0)]
  477.     pln = crossp(p0, pole)
  478.     split=[!map.p0lon, !map.p0lat, pln, 0.0]
  479. endelse
  480.  
  481. noadd = keyword_set(add) eq 0    ;= 0 if we're adding
  482. MAP_CLIP_SET, RESET=noadd, SPLIT=split, TRANSFORM=noadd
  483. end
  484.  
  485.  
  486.  
  487. PRO map_set_clip, pname
  488. ; Set up the default clipping/splitting for the given projection.
  489. ;
  490. r = -1.0            ;Great circle clipping radius...
  491. if pname eq "LambertConic" or pname eq "AlbersEqualAreaConic" then begin
  492.     isign = 2*(!map.p[3] ge 0.0)-1 ;Assume its in one hemisphere
  493. ;    map_set_split, /EQUATOR
  494.     map_set_split
  495.     map_clip_set, CLIP_PLANE=[0,0,isign,0]
  496.     map_clip_set, CLIP_PLANE=[0,0,- isign, .95] ;Towards pole
  497.     
  498. endif else if pname eq 'Mercator' or pname eq 'MillerCylindrical' then begin
  499.     rm = 0.975                  ;Go out 97.5% of the way to the poles
  500.     map_set_split
  501.     map_clip_set, CLIP_PLANE=[-!map.pole[4:6], rm]
  502.     map_clip_set, CLIP_PLANE=[!map.pole[4:6], rm]
  503.     
  504. endif $                         ;Other cylindricals
  505. else if pname eq "Stereographic" or pname eq "Orthographic" then r = -1.0e-5 $
  506. else if pname eq "Satellite" then r = - 1.01/!map.p[0] $
  507. else if pname eq "Gnomic" then r = -0.5 $
  508. else if pname eq "LambertAzimuthal" or pname eq "AzimuthalEquidistant" then $
  509.   r=0.8 $
  510. else if pname eq "TransverseMercator" then r = -0.4 $ ;UTM, a fudge factor.
  511. else begin                      ;If we haven't setup a clip, add splitting
  512.     map_proj_info, /CURRENT, CYLINDRICAL=is_cyl
  513.     if is_cyl then map_set_split
  514. endelse
  515.  
  516. ;    If r is set, set clipping to points on a plane whose normal
  517. ;    passes thru the center (u0,v0), and at a distance of r from the
  518. ;    origin.  r < 0 is the side towards (u0,v0).
  519.  
  520. if r ne -1.0 then $
  521.   map_clip_set, CLIP_PLANE=[cos(!map.u0) * !map.coso, $
  522.                             sin(!map.u0) * !map.coso, !map.sino, r]
  523. map_clip_set, /TRANSFORM
  524. end
  525.  
  526.  
  527.  
  528. ; *****************************************************************
  529. PRO MAP_SET, p0lat, p0lon, rot, $
  530. ;    ********** Projection keywords:
  531.   PROJECTION=proj, $              ;The projection index
  532.   NAME=name, $                    ;The projection name as a string
  533.   STEREOGRAPHIC = stereographic, $
  534.   ORTHOGRAPHIC = orthographic, $
  535.   CONIC = conic, $
  536.   LAMBERT = lambert, $
  537.   GNOMIC = gnomic, $
  538.   AZIMUTHAL = azimuth, $
  539.   SATELLITE = satellite, $      ;Also called general perspective proj
  540.   CYLINDRICAL = cylindrical, $
  541.   MERCATOR = mercator, $
  542.   MILLER_CYLINDRICAL=miller, $
  543.   MOLLWEIDE = mollweide, $     
  544.   SINUSOIDAL = sinusoidal, $
  545.   AITOFF = aitoff, $            ;Original Aitoff projection
  546.   HAMMER = hammer, $            ;This is really the Hammer-Aitoff projection
  547.   ALBERS = albers, $            ;Alber's equal-area conic
  548.   TRANSVERSE_MERCATOR = utm, $  ;Transverse mercator for an ellipsoid
  549.   ROBINSON = robinson, $        ;Robinson projection
  550.                                 ; Also called Gauss-Krueger in Europe.
  551. ; **** Projection specific keywords:
  552. ELLIPSOID = ellips, $           ;Defines the ellipsoid for Transverse mercator.
  553.                                 ;3 elements: a = equatorial radius (meters), e^2
  554.                                 ;= eccentricity squared, k0 = scale on
  555.                                 ;central meridian. e^2 = 2*f-f^2,
  556.                                 ;where f = 1-b/a, b = polar radius
  557.                                 ;Default is the Clarke 1866 ellipsoid, =
  558.                                 ; [6378206.4d0, 0.00676866d0, 0.9996d0]
  559. CENTRAL_AZIMUTH=cent_azim, $    ;Angle of central azimuth (degrees) for:
  560.                                 ; Cylindrical, Mercator, Miller,
  561.                                 ; Mollw, Robinson, and Sinusoidal
  562.                                 ; projections. Default = 0.
  563.                                 ; The pole is placed at an azimuth of
  564.                                 ; CENTRAL_AZIMUTH degrees CCW of north.
  565. STANDARD_PARALLELS = std_p, $   ;One or two standard parallels for conics,
  566.                                 ; One or two element array.
  567. SAT_P = Sat_p, $                ;Satellite parameters: Altitude expressed in
  568.                                 ; units of radii, Omega, and rotation.
  569.                                 ; Rotation may also specified by the
  570.                                 ; rot parameter.
  571. ;    ********** MAP_SET specific keywords:
  572. CLIP=clip, $                    ;Default = do map specific clipping,
  573.                                 ; CLIP=0 to disable
  574. SCALE=scale, $                  ;Construct isotropic map with a given scale.
  575.                                 ; Map Scale is 1:scale, otherwise fit to window
  576. ISOTROPIC = iso, $,             ;If set, make X and Y scales equal
  577. LIMIT = limit, $                ;4 or 8 point lat/lon limit:
  578.                                 ; 4 point: [latmin, lonmin, latmax, lonmax]
  579.                                 ; 8 point: [latLeft,lonLeft, latTop,
  580.                                 ;       lonTop, LatRt, lonRt, LatBot, LonBot]
  581. ;    ********** MAP_SET graphics keywords:
  582. NOERASE=noerase, TITLE=title,$
  583.   ADVANCE = advance, COLOR=color, POSITION = position, $
  584.   NOBORDER=noborder, T3D=t3d, ZVALUE=zvalue, $
  585.   CHARSIZE = charsize, XMARGIN=xmargin, YMARGIN=ymargin, $
  586. ;    ********** MAP_HORIZON keywords:
  587. HORIZON=horizon, E_HORIZON=ehorizon, $ ; E_HORIZON = structure containing
  588.                                 ; extra keywords passed to the
  589.                                 ; map_horizon procedure, e.g.
  590.                                 ; E_HORIZON={FILL:1}
  591. ;    ********** MAP_CONTINENTS keywords:
  592. CONTINENTS = continents, E_CONTINENTS=econt, $ ;E_CONTINENTS = structure
  593.                                 ; containing extra keywords to be
  594.                                 ; passed to
  595.                                 ; map_continents, e.g. E_CONTINENTS={FILL:1}
  596.   USA=usa, HIRES = hires, $
  597.   MLINESTYLE=mlinestyle, MLINETHICK=mlinethick, CON_COLOR=con_color, $
  598.  
  599. ;    ********** MAP_GRID keywords:
  600.   GRID=grid, E_GRID=egrid, $    ;E_GRID = extra keywords structure
  601.   GLINESTYLE=glinestyle, GLINETHICK=glinethick, $
  602.   LABEL=label, LATALIGN=latalign, LATDEL=latdel, LATLAB=latlab, $
  603.   LONALIGN=lonalign, LONDEL=londel, LONLAB=lonlab, $
  604.  
  605. ; Ignored, but here for compatibility:
  606.   WHOLE_MAP=whole_map
  607. ; *****************************************************************
  608. ;+
  609. ;         Limit         =  A four or eight element vector.
  610. ;              If a four element vector, [Latmin, Lonmin, Latmax,
  611. ;                          Lonmax] specifying the boundaries of the region
  612. ;                          to be mapped. (Latmin, Lonmin) and (Latmax,
  613. ;                          Lonmax) are the latitudes and longitudes of
  614. ;                          two diagonal points on the boundary with
  615. ;                          Latmin < Latmax and Lonmin < Lonmax.
  616. ;               For maps that cross the international dateline,
  617. ;               specify west longitudes as angles between
  618. ;               180 and 360.
  619. ;              If an eight element vector: [ lat0, lon0, 
  620. ;                lat1, lon1, lat2, lon2, lat3, lon3] specify
  621. ;                four points on the map which give, respectively,
  622. ;                the location of a point on the left edge, 
  623. ;                top edge, right edge, and bottom edge of
  624. ;                the map extent.
  625. ;-
  626.  
  627.  
  628. ; If the map range can be simply expressed as [latmin, lonmin, latmax,
  629. ; lonmax], then they are contained here.  With many maps, this is
  630. ; impossible. These values are used to speed MAP_GRID and
  631. ; MAP_CONTINENTS, by not drawing elements that are clearly off the
  632. ; map. 
  633. common map_set_common
  634.  
  635. on_error, 2                     ;Return to caller if error.
  636. del = 1.0e-6
  637. dtor = !dpi/ 180.
  638. dpi2 = !dpi/2.
  639. !map.ll_box = 0                 ;Assume no limits
  640.  
  641. if n_elements(proj) then iproj = proj $
  642. else if n_elements(name) then begin ;Name supplied, match it
  643.     map_proj_info, PROJ_NAMES = names ;Available names
  644.     name1 = strcompress(strupcase(name), /REMOVE_ALL) ;Ignore case & blanks
  645.     iproj = (where(name1 eq strupcase(names)))(0) ;Match it
  646.     if iproj lt 1 then message, 'Projection "'+name+'" not available'
  647. endif else if keyword_set(stereographic) then iproj =1 $
  648.   else if keyword_set(orthographic) then iproj =2 $
  649.   else if keyword_set(conic  ) then iproj =3 $
  650.   else if keyword_set(lambert) then iproj =4 $
  651.   else if keyword_set(gnomic ) then iproj =5 $
  652.   else if keyword_set(azimuth) then iproj =6 $
  653.   else if keyword_set(satellite) then iproj =7 $
  654.   else if keyword_set(mercator) then iproj =9 $
  655.   else if keyword_set(mollweide) then iproj = 10 $
  656.   else if keyword_set(sinusoidal) then iproj = 11 $
  657.   else if keyword_set(aitoff) then iproj = 12 $
  658.   else if keyword_set(hammer) then iproj = 13 $
  659.   else if keyword_set(albers) then iproj = 14 $
  660.   else if keyword_set(utm) then iproj = 15 $
  661.   else if keyword_set(miller) then iproj = 16 $
  662.   else if keyword_set(robinson) then iproj = 17 $
  663.   else iproj=8                                     ;Assume cylindrical
  664.  
  665. !map.projection = iproj
  666. map_proj_info, /CURRENT, NAME=pname ;Get symbolic name
  667.  
  668. ; Initial erase?
  669. if !P.multi[0] eq 0 and keyword_set(advance) then erase
  670. if (not KEYWORD_SET(noerase)) and (not KEYWORD_SET(advance)) THEN erase
  671.  
  672.  ; ******** Supply defaults  **********
  673.  
  674. if n_params() lt 3 then rot = 0.0d0
  675. if n_params() lt 2 then p0lon = 0d0
  676. if n_params() lt 1 then p0lat = 0d0
  677.  
  678. if n_elements(cent_azim) le 0 then cent_azim = 0.0
  679. if n_elements(color) eq 0 then color = !p.color  ;Default color
  680. if n_elements(title) eq 0 then title = " "
  681. if n_elements(t3d) le 0 then t3d = 0
  682. if n_elements(zvalue) eq 0 then zvalue = 0.
  683. if n_elements(charsize) eq 0 then charsize = !p.charsize
  684. if charsize le 0.0 then charsize = 1.0
  685. if n_elements(clip) eq 0 then clip = 1 ;default for CLIP is ON
  686.  
  687.                                 ;Hack for CONIC backwards compatibility
  688. if pname eq "LambertConic" and (rot ne 0) and n_elements(std_p) eq 0 then begin
  689.     std_p = [p0lat, p0lon]      ;Two standard parallels
  690.     p0lon1 = rot                ;and rotation is the longitude
  691. endif else p0lon1 = p0lon
  692.  
  693. if abs(p0lat) gt 90.0 then $
  694.    message,'Latitude must be in range of +/- 90 degrees'
  695. if abs(p0lon) gt 360.0 then $
  696.    message,'Longitude must be in range of +/- 360 degrees'
  697.  
  698. u = double(p0lon1)               ;Reduce lon to +/- 180
  699. while u lt -180.0 do u = u + 360.
  700. while u gt 180.0 do u = u - 360.
  701.  
  702. !map.p0lon = u
  703. !map.p0lat = p0lat
  704. !map.u0 = u * dtor
  705. !map.v0 = p0lat * dtor
  706.  
  707. !x.type = 3                     ;New map code
  708. !y.type = 0
  709.  
  710. ;    Simple projection?  (used for (psuedo) cylindricals)
  711. !map.simple = abs(p0lat) le del and abs(cent_azim) le del
  712.  
  713. !map.sino = sin(!map.v0)
  714. !map.coso = cos(!map.v0)
  715.  
  716. !map.rotation = rot
  717. !map.sinr = sin(rot * dtor)
  718. !map.cosr = cos(rot * dtor)
  719.  
  720. ;    Compute position of Pole which is a distance of !PI/2 from
  721. ;    (p0lon, p0lat),    at an azimuth of rot CCW of north.
  722. ;
  723. pole = [0.0, sin(cent_azim * dtor), cos(cent_azim * dtor)]
  724. p2 = map_rotxyz(pole, 0, -p0lat, u)    ;Rotate to put origin at (0,1,0)
  725. plat = asin(p2[2])
  726. cosla = sqrt(1.0-p2[2]^2)
  727. if cosla eq 0.0 then begin         ;On pole?
  728.     plon = 0.0 & sinln = 0.0 & cosln = 1.0
  729.   endif else begin
  730.     plon = atan(p2[1], p2[0])
  731.     sinln = p2[1] / cosla & cosln = p2[0] / cosla
  732.   endelse
  733.  
  734.         ; lon/lat, sin(lat), cos(lat), xyz, Location of pole
  735. !map.pole = [plon, plat, p2[2], cosla, p2]
  736.  
  737. if pname eq "LambertConic" or pname eq "AlbersEqualAreaConic" then begin
  738.     if n_elements(std_p) eq 0 then begin
  739.         if p0lat ne 0 then s = [p0lat, p0lat] else s = [20., 20.]
  740.     endif else if n_elements(std_p) eq 1 then s = std_p[[0,0]] $
  741.     else if n_elements(std_p) eq 2 then s = std_p $
  742.     else message, 'STANDARD_PARALLELS must have 1 or 2 elements'
  743.     s = s * dtor
  744.     if pname eq "LambertConic" then begin
  745.         if ABS(!map.p0lat) eq 90. then $
  746.           message, 'Center of projection may not be pole for conics'
  747.         if s[0] * s[1] lt 0 then $
  748.           message, 'Standard parallels must be in same hemisphere'
  749.         if s[0] eq s[1] then n = sin(s[0]) $
  750.         else n = alog(cos(s[0])/cos(s[1])) / $
  751.           alog(tan(!dpi/4+s[1]/2)/tan(!dpi/4+s[0]/2))
  752.         if n eq 0.0 then message, $
  753.           'Equator is illegal standard parallel for conic'
  754.         F = cos(s[0]) * tan(!dpi/4 + s[0]/2)^n/n
  755.         r0 = F/tan(!dpi/4+!map.v0/2.)^n
  756.         !map.p = [n, F, r0, s]
  757.     endif else begin            ;Albers conic
  758.         n = (sin(s[0]) + sin(s[1]))/2.
  759.         if n eq 0 then begin
  760.             message,'Standard parallels should not be 0', /INFO
  761.             n = 0.5
  762.         endif
  763.         c = cos(s[0])^2 + 2 * n * sin(s[0])
  764.         r0 = sqrt(c - 2*n*sin(!map.v0))/n
  765.         !map.p = [n, c, r0, s]
  766.     endelse
  767. endif
  768.  
  769. if pname eq "Satellite" then begin ;Special params for satellite projection.
  770.     if n_elements(sat_p) eq 0 then !map.p[0] = 2.0 $ ;Altitude in radii
  771.     else !map.p[0] = sat_p[0]
  772.     if n_elements(sat_p) le 1 then omega = 0d0 else omega = dtor * sat_p[1]
  773.                                 ;Save em.  sat(1) = TRUE for simple case (Vertical perspective)
  774.     !map.p[1] = omega
  775.     !map.p[2] = sin(omega)    ;Somega = p[1]
  776.     !map.p[3] = cos(omega)    ;comega
  777.     if n_elements(sat_p) ge 3 then begin ; For backwards compatibility, if
  778. ; Sat_p contains three elements, interpret the 3rd element as the rotation:
  779.         !map.rotation = sat_p[2]
  780.         !map.sinr = sin(sat_p[2] * dtor)
  781.         !map.cosr = cos(sat_p[2] * dtor)
  782.     endif                       ;3 elements
  783. endif
  784.  
  785.  
  786. map_proj_info, /CURRENT, WIDTH=meters, UVRANGE=uvrange, $
  787.   CYLINDRICAL = is_cyl
  788.  
  789. ; Set up cylindrical projections
  790. if pname eq "TransverseMercator" then begin ;Special params for UTM
  791.     if n_elements(ellips) ne 3 then $ ;Default == Clarke 1866 ellipsoid
  792.       ellips = [6378206.4d0, 0.00676866d0, 0.9996d0]
  793. ;   ellips is a 3 element array containing the ellipsoid parameters:
  794. ;           [a, e^2, k0]
  795. ;    a = Equatorial radius, in meters.
  796. ;    b = polar radius, b = a * (1-f), f = 1-b/a
  797. ;    e^2 = eccentricity^2 = 2*f-f^2, where f = flattening.
  798. ;    k0 = scale on central meridian, = 0.9996 for UTM.
  799.     e_2 = ellips[1]/(1.0-ellips[1])
  800. ;        k0      e2       e_2   a         b      m0  m
  801.     !map.p = [ellips[2], ellips[1], e_2, ellips[0], ellips[1], 0., 0.]
  802.     !x.s = [0,1] & !y.s = [0,1] & !x.type = 3 ;Fake coordinate system.
  803.     q = convert_coord(u, p0lat, /data, /to_norm) ;Get m0
  804.     !map.p[5] = !map.p[6]     ;Set m0 from Center of projection
  805.     
  806. endif else if is_cyl then begin ;other Cylindrical or pseudo-cylind proj?
  807. ; I don't know why we can't solve for this angle (the azimuth of (u0,v0)
  808. ; from (xp, yp)) using the law of sines.  This is solving it the hard
  809. ; and long way.....  But it works....
  810.     az = atan(!map.coso * sin(!map.u0-plon), $
  811.               cosla * !map.sino - p2[2] * !map.coso  *cos(!map.u0-plon))
  812.     !map.p[0] = dpi2-az
  813. ;    print, 'Pole: ', plon/dtor, plat /dtor, ', Az=',az / dtor
  814. endif
  815.  
  816.  
  817. uvrange_orig = uvrange          ;Save original UV
  818.  
  819. if is_cyl then begin            ;Adjust for rotation & wrap in cylind projs
  820. ;    Set the x coordinate of the edge in UV space for wrap-around detection.
  821. ;    This is only revelant in (pseudo) cylindrical projections.
  822.     u = uvrange[[0,0,2,2]]    ;The 4 corners.
  823.     v = uvrange[[1,3,1,3]]
  824.     t1 = u * !map.cosr + v * !map.sinr ;Now rotate by - angle...
  825.     v =  v * !map.cosr - u * !map.sinr
  826.     uvrange[0] = min(t1, max=t2)
  827.     uvrange[2] = t2
  828.     uvrange[1] = min(v,  max=t2)
  829.     uvrange[3] = t2
  830. endif
  831.  
  832. if n_elements(limit) gt 0 then begin ;Get UV range.
  833.     if keyword_set(scale) then $
  834.       message, 'Conflicting keywords specified: LIMIT and SCALE', /INFO
  835.     lim = limit                 ;Save old
  836.     if n_elements(lim) eq 4 then !map.ll_box = lim ;Save if 4 point limit
  837.     map_set_limits, lim, uvrange
  838. endif
  839. ; print,'Uvrange=', uvrange
  840.  
  841. ;*****************************************************************
  842.  
  843. if keyword_set(position) then  $
  844.     plot, [0,1], xsty=5, ysty=5, /nodata,    $
  845.     title=title, /noerase, position=position, charsize=charsize $
  846. else begin            ;Position not set
  847.                                 ; Use margins
  848.     if n_elements(xmargin) eq 1 then xmar = [xmargin, xmargin] $
  849.     else if n_elements(xmargin) eq 2 then xmar = xmargin $
  850.     else xmar = [1,1]
  851.     if n_elements(ymargin) eq 1 then ymar = [ymargin, ymargin] $
  852.     else if n_elements(ymargin) eq 2 then ymar = ymargin $
  853.     else ymar = [1,2]
  854.     plot, [0,1], XSTYLE=5, YSTYLE=5, /NODATA, COLOR=color,    $
  855.       XMARGIN= xmar, YMARGIN= ymar, TITLE=title, /NOERASE, $
  856.       CHARSIZE = charsize
  857. endelse
  858.  
  859. !x.type = 3            ;For new maps;
  860.  
  861. if keyword_set(noborder) eq 0 then begin
  862. ; fudge is used to add a bit of spacing between the border and the
  863. ; extent of the map region. fudge of 0.01 indicates that the extra
  864. ; spacing (internal map margin) should be 1% of the original map extent.
  865. ;
  866. ; this extra space is now turned off when the noborder keyword is set
  867. ;
  868.     fudge = 0.01
  869.     s = (uvrange[2]-uvrange[0]) * fudge
  870.     ss = (uvrange[3] - uvrange[1]) * fudge
  871.     uvrange = uvrange + [-s, -ss, s, ss]
  872. endif
  873.  
  874. ;     Figure the size of the drawing area on the screen in normalized units
  875. x_size = !x.window[1]-!x.window[0]
  876. y_size = !y.window[1]-!y.window[0]
  877. x_cen = total(!x.window)/2.     ;Midpoints in X & Y
  878. y_cen = total(!y.window)/2.
  879.  
  880. ; Compute the X and Y scale factors, !x.s(1), !y.s(1):
  881. if keyword_set(scale) then begin ;Absolute scale provided?
  882.                                 ;Absolute size of entire map, meters / uvrange
  883.     if meters ne 0 then $       ;Anything with an inherent scale 
  884.       s = (meters / scale) / (uvrange_orig[2]-uvrange_orig[0])$ ;meters/uv_unit
  885.     else s = 1.0 / scale        ;projections already in meters
  886.     
  887.     ss = !d.x_size / !d.x_px_cm / 100. ;width of window in meters
  888.     !x.s[1] = s / ss
  889.     !y.s[1] = !x.s[1] * !d.x_size / !d.y_size ;Adjust for aspect ratio
  890.             ;Adjust the UV range to include the plotting window  **** 
  891.     t = [x_size / !x.s[1], y_size / !y.s[1]] / 2 ;half the UV range
  892.     c = [uvrange[0]+uvrange[2], uvrange[1]+uvrange[3]]/2 ;UV center
  893. ; Take smallest UV box:
  894.     uvrange = [c[0]-t[0] > uvrange[0], c[1]-t[1] > uvrange[1], $
  895.                c[0]+t[0] < uvrange[2], c[1]+t[1] < uvrange[3]] ;New uv range
  896.  
  897. endif else if keyword_set(iso) then begin
  898.  ; Scale in pixels/uvunit, correct for ISOMETRIC aspect ratio. Use smaller
  899.  ; of X and Y scale factors
  900.     sx = x_size * !d.x_size /(uvrange[2]-uvrange[0]) ;X scale
  901.     sy = y_size * !d.y_size/(uvrange[3]-uvrange[1]) ;Y scale
  902.  
  903.     !x.s[1] = (sx < sy) / !d.x_size ;Set smaller
  904.     !y.s[1] = !x.s[1] * !d.x_size / !d.y_size
  905.     if sx gt sy then $          ;Resize window to fit map area
  906.       !x.window = x_cen + sy/sx/2 * [-x_size, x_size] $
  907.     else !y.window = y_cen + sx/sy/2 * [-y_size, y_size]
  908. endif else begin                ;Scale to fit WINDOW
  909.     !x.s[1] = x_size/(uvrange[2]-uvrange[0])
  910.     !y.s[1] = y_size/(uvrange[3]-uvrange[1])
  911. endelse
  912.  
  913. ; Compute offsets to center UV rectangle in the window
  914. !x.s[0] = -(uvrange[0]+uvrange[2])/2. * !x.s[1] + x_cen
  915. !y.s[0] = -(uvrange[1]+uvrange[3])/2. * !y.s[1] + y_cen
  916.  
  917. map_clip_set,/reset        ;Clear clipping pipeline.
  918.  
  919. if keyword_set(clip) then begin    ;Setup clipping
  920.     map_set_clip, pname         ;Set default clipping?
  921.     if (n_elements(lim) gt 0) or keyword_set(scale) then $ ;Clip UV?
  922.       MAP_CLIP_SET, CLIP_UV = uvrange
  923.                                 ; Set up plotting clip window
  924.     !p.clip([0,2]) = (uvrange([0,2]) * !x.s(1) + !x.s(0)) * !d.x_size
  925.     !p.clip([1,3]) = (uvrange([1,3]) * !y.s(1) + !y.s(0)) * !d.y_size
  926. endif
  927.  
  928. !map.uv_box = uvrange          ;Save UV range
  929.  
  930. if clip and total(abs(!map.ll_box)) eq 0 then $ ; Try to make a lon/lat range
  931.   map_set_ll_box
  932.  
  933. if keyword_set(noborder) eq 0 then  $ ;Draw the border.
  934.   plots, !x.window[[0,1,1,0,0]], !y.window[[0,0,1,1,0]], $
  935.      COLOR=color, zvalue, /NORM, /NOCLIP, T3D=t3d
  936.  
  937. ; Collect the common graphics keywords
  938. map_struct_append, egraphics, "COLOR", color
  939. map_struct_append, egraphics, "T3D", t3d
  940. map_struct_append, egraphics, "ZVALUE", zvalue
  941.  
  942. ; Process MAP_HORIZON keywords:    **************
  943. map_struct_merge, ehorizon, egraphics ;Add common graphics keywords
  944. if keyword_set(horizon) then MAP_HORIZON, _EXTRA=ehorizon
  945.  
  946. ; Process MAP_CONTINENT keywords:    **************
  947. if n_elements(mlinestyle) then map_struct_append, econt, "LINESTYLE", mlinestyle
  948. if n_elements(mlinethick) then map_struct_append, econt, "THICK", mlinethick
  949. if n_elements(con_color) then map_struct_append, econt, "COLOR", con_color
  950. if n_elements(hires) then map_struct_append, econt, "HIRES", hires
  951. if n_elements(continents) then map_struct_append, econt, "CONTINENTS", continents
  952. if n_elements(usa) then map_struct_append, econt, "USA", usa
  953. if n_elements(econt) gt 0 or keyword_set(continents) then begin
  954.     map_struct_merge, econt, egraphics ;Add common graphics kwrds
  955.     MAP_CONTINENTS, _EXTRA=econt
  956.     endif
  957.  
  958. ; Process MAP_GRID keywords:    **************
  959. if n_elements(label)  then map_struct_append, egrid, "LABEL", label
  960. if n_elements(latlab) then map_struct_append, egrid, "LATLAB", latlab
  961. if n_elements(lonlab) then map_struct_append, egrid, "LONLAB", lonlab
  962. if n_elements(latdel) then map_struct_append, egrid, "LATDEL", latdel
  963. if n_elements(londel) then map_struct_append, egrid, "LONDEL", londel
  964. if n_elements(latalign) then map_struct_append, egrid, "LATALIGN", latalign
  965. if n_elements(lonalign) then map_struct_append, egrid, "LONALIGN", lonalign
  966. do_grid = (keyword_set(grid) + n_elements(egrid) + n_elements(glinestyle) + $
  967.            n_elements(glinethick)) ne 0
  968.  
  969. map_struct_merge, egrid, egraphics
  970. if n_elements(glinestyle) then $ ;These keywords supersede those in egraphics
  971.   map_struct_append, egrid, "LINESTYLE", glinestyle
  972. if n_elements(glinethick) then $
  973.   map_struct_append, egrid, "THICK", glinethick
  974. if do_grid then MAP_GRID, CHARSIZE=charsize, _EXTRA=egrid
  975.  
  976. if KEYWORD_SET(ADVANCE) and !P.Multi[0] gt 0 THEN $
  977.      !P.Multi[0] = !P.Multi[0] - 1 $
  978. else !p.multi[0] = !p.multi[1] * !p.multi[2] - 1
  979.  
  980. end
  981.